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ABSTRACT 

We present the first results of a global axisymmetric simulation of accretion onto ro- 
tating magnetized stars from a turbulent accretion disk, where the turbulence is driven 
by the magnetorotational instability (MRI). Long-lasting accretion is observed for sev- 
eral thousand rotation periods of the inner part of the disc. The angular momentum 
is transported outward by the magnetic stress of the turbulent flow with a rate corre- 
sponding to a Shakura-Sunyaev viscosity parameter a ^ 0.01 - 0.04. Close to the star 
the disk is stopped by the magnetic pressure of the magnetosphere. The subsequent 
evolution depends on the orientation of the poloidal field in the disk relative to that of 
the star at the disk-magnetosphere boundary. If fields have the same polarity then the 
magnetic flux is accumulated at the boundary and blocks the accretion which leads 
to the accumulation of matter at the boundary. Subsequently this matter accretes to 
the star in outburst before accumulating again. Hence, the cycling, 'bursty' accretion 
is observed. The magnetic stress is enhanced at the boundary leading to the enhanced 
accretion rate. If the disc and stellar fields have opposite polarity then the field re- 
connection enhances the penetration of the disk matter towards the deeper field lines 
of the magnetosphere. However, the magnetic stress at the boundary is lower due to 
the field reconnection. This decreases the accretion rate and leads to smoother accre- 
tion at a lower rate. Test simulations show that in the case of higher accretion rate 
corresponding to a = 0.05 — 0.1, accretion is bursty in cases of both polarities. On 
the other hand, at much lower accretion rates corresponding to a < 0.01, accretion 
is not bursty in any of these cases. We conclude that the episodic, bursty accretion is 
expected during periods of higher accretion rates in the disc, and in some cases it may 
alternate between bursty and smooth accretion, if the disk brings the poloidal field 
of alternating polarity. We find that a rotating, magnetically-dominated corona forms 
above and below the disk, and that it slowly expands outward, driven by the magnetic 
force. 
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1 INTRODUCTION 

The magnetorotational instability (MRI) is the likely origin 
of turbulent stress in accretion disks around black holes and 
stars (e.g., [Balbus & Hawley||1991| [Balbus & Hawley||1998| ). 
MRI-driven accretion has been observed in a number of lo- 
cal and global axisymmetric and 3D MHD simulations (e.g., 



Hawley et al. 1995 Stone et al. 1996; Gammie & Menou 



1998, ,Stone & Pringle||2001[ |Hawley..2QQQ, .Hawley et al. 
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|200H [Hawley & Krolik|2 002VBec kwith, et al.|20 09). In these 
simulations the central object is a black hole. However, there 
are many stars which have dynamically important magnetic 
fields. These include young, Classical T Tauri stars (hereafter 
CTTSs; e.g., |Bouvier et al.|2005| ), some types of white dwarfs, 
and neutron stars in binary systems (e.g., Warner|2004[ |Van| 
[der Klis 2000 ). In these stars, the magnetic field is strong 
enough to stop the disk at radii larger than the radius of 
the star. Many observational properties are determined by the 
processes at the disk-magnetosphere boundary. 

The disk-magnetosphere interaction in the case of lam- 
inar (non-turbulent) flow in the disk has been investigated 
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Figure 1. Initial distribution of the magnetic field used in different models. The dipole magnetic field of the star is fixed at Jl = 20, but the 
field distribution in the disc is different. From left to right: the field in the disc is vertical (parallel to the z— axis) and has the same direction 
as the magnetic field of the star at the disc-magnetosphere boundary (model VP20); same but for the opposite direction of the field in the disc 
(antiparallel case, model VA20); tapered field with z— component parallel to the field of the star (model TP20); tapered field with antiparallel 
field (model TA20); loops in the disc where the inner loop field is parallel to the star's field (model LP20). 




Figure 2. Left panel: The grid used in simulations; Middle panel: same 
grid, but rarefied by a factor of 10. Right panel: Sample simulation run 
shows the density distribution and selected field lines at the end of a 
long simulation run in model TP20 att = 3000. 
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|Long et al. 2068\. The laminar flow in the 



|& Sunyaev|1973p . Similarly, the a— type magnetic diffusivity 
has been included in a number of studies. Heretofore, accre- 
tion due to turbulent, MRI-driven discs has not been studied. 

Here, we present the results of the first axisymmetric 
simulations of MRI-driven accretion onto magnetized stars. 
In these simulations we were able to use a very high grid 
resolution and perform very long runs. Results for global 3D 
simulations will be discussed in a separate paper ( [Romanova] 
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Table 1. The Table describes the main simulation models, where the 
first letter in the model name stands for the field configuration: 'V - 
vertical field, 'T' - tapered field, 'L' - looped field. The second letter 
stands for the mutual directions of the magnetic fields of the disc and 
the star: 'P' - parallel and A' - antiparallel. The number at the end of 
the name shows the value of the dimensionless magnetic moment ji. 
The Table also shows the z— component of the seed magnetic field in 
the disc, B^, and the duration of simulations in periods of rotation 
Pmid at the grid center, r = 22 (Pmid = IOSPq)- 



|et al.|20lT ) . Accretion onto rapidly rotating stars will be de- 
scribed in Ustyugova et al.| ( |2011| ). Here, we consider the case 
of slowly rotating stars. 
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2 THEORETICAL BACKGROUND 

The magnetorotational instability (MRI) arises under con- 
ditions where a weak magnetic field threads the disc. Be- 
low, we briefly summarize the condition for the onset of the 
MRI instability (Balbus & Hawley 1991) for a simple case 
where an axial magnetic field Bqz threads a Keplerian disc 
which rotates with angular rate Q = {GM/r^y^'^ where 
M is the mass of the central object. For axisymmetric per- 
turbations of the disc with 5v = [6vr{z,t),6vcf>{z,t),0] and 
= [6Br{z,t),6Bci){z,t),0] and for perturbations propor- 
tional to e:>^p{ikzZ — iut), one finds the dispersion relation 
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where 



va = Bq/ y^4np is the Alfven velocity and Kr = [41^^ + 
2rQdQ/dr]-^^'^ is the radial epicycHc frequency of the disc. In 
order for the perturbation to fit within the vertical extent of 
the disc one needs kzh > 1, where h = Cg/Q is the half- 
thickness of the disc and Cs is the isothermal sound speed in 
the disc. For most conditions the disc is thin with h <^ r, or 
Cs < rQ. 

Evidently there can be instability if cj?. < which hap- 
pens if (kvA)'^ < —2rQdQ/dr. For a Keplerian disc this cor- 
responds to (kvA)'^ < 3Q^. Therefore, the above-mentioned 
condition that kzh> 1 implies that the instability occurs only 

for VA < Cs, or 



2ci 



> 1 



(2) 



Note that /3 is based on the initial vertical magnetic field, 
Bq. As a result of the instability the magnetic field may 
grow to values much larger than Bq. The maximum value 
of the growth rate is ^(cj)max = 3Q/4, and it occurs for 
A^max = (15/16)^/^Q/'yA. For /3 < 1 the perturbation does 
not fit inside the disc and there is stability. As /3 increases from 
unity (weaker magnetic field) the maximum growth rate stays 
the same but the wavelength of the perturbation gets shorter 
(oc /3~^/^). Of course, for sufficiently small wavelengths the 
damping due to numerical viscosity z^num/c^) will be larger 
than the MRI growth rate. To generate MRI-driven accretion 
in numerical simulations, one should have large enough /3 in 
the disc so as to have instability, and at the same time high 
enough grid resolution to avoid the damping of small wave- 
lengths. 



3 SIMULATIONS OF MRI INSTABILITY 
3.1 The problem setup 

We investigate MRI-driven accretion onto a rotating star with 
a dipole magnetic field. The field lines are frozen into the 
stellar surface, and initially, they thread the entire simulation 
region, including the disc. The disc is dense and cold while 
the corona is 10^ times hotter and 10^ times less dense. The 
disc rotates with the Keplerian velocity at the beginning of 
the simulation. The corona also rotates with Keplerian veloc- 
ities at cylindrical radii corresponding to the rotation of the 
disc. This helps to avoid the build-up of an initial discontinu- 
ity of the magnetic field at the disc-magnetosphere boundary. 
We derive the initial distribution of density and pressure in 



the disc and corona by balancing the gravitational, centrifu- 
gal and pressure gradient forces (see Sec.|A4|for details). The 
initial inner disc radius is set to r = 10 (in units of stellar 
radius, Ro; see Sec. |m] for the full description of reference 
units) . The star rotates slowly with an angular velocity corre- 
sponding to a corotation radius (rcor = (GM/Qiy^^) of 10 
(in dimensionless units) . 

A small seed poloidal field with 2;— component Bd has 
been added to the disc, which is necessary for the initial ex- 
citation of the MRI-driven turbulence. We considered three 
types of initial configurations: 

• Vertical fieldiV-' type models): The field is vertical (that is, 
parallel to the 2;— axis), has a value of Bd and is homogeneous 
in space. The alignment of the field may be parallel or antipar- 
allel to the star's magnetic field at the disc-magnetosphere 
boundary (see two left panels in Fig.[T]). 

• Tapered field ('T-'type models) : An initially vertical field 
tapered in such a way that the magnetic flux is restricted to 
lie within the disc . The magnetic flux inside a disc with half- 
thickness /i (at 1^1 < h) is determined by: 
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where ^c(r) — K.GM/r, E is the constant of integration in 
the initial equilibrium equation (see Sec. |A4| ). Fig. [T] (3rd and 
4th panels from the left) shows cases of parallel (model 'TP') 
and antiparallel (model 'TA') fields. 

• Looped field ('L-'type models): The magnetic flux inside a 
disc with half-thickness h is determined by: 



Bdv' ( z 
— - — cos avr — - 
2 \ 2h 
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■ Td 



Tout - Td 



(4) 



where parameters a and h determine the number of loops in 
the vertical and horizontal directions, rout is the radius of the 
external boundary, and is the inner disc radius. 

We experimented with different numbers of loops with val- 
ues of A^rioops X A^^ioops ranging from 3 x 1 to 7 x 5, and con- 
cluded that the case with the smaller number of loops (3x1) 
is sufficient for exciting the turbulence, while the cases with 
a larger number of loops do not have advantages. The total 
poloidal flux is not zero, and hence we have two cases where 
the averaged field in the disc is parallel (model 'LP') or an- 
tiparallel (model 'LA') to the field of the star (see an example 
for 'LP-'type configuration at the right panel of Fig. [T]) . 



We solve the ideal MHD equations (see Sec. |A1| ), given 
axisymmetric conditions, in cylindrical coordinates using a 
Godunov-type code (see Sec. |A3| ). The grid is of high reso- 
lution with compression towards the disc midplane and to- 
wards the axis with the goal of having a larger number of 
the grid-points in the disc and near the star (see Fig. [2| . After 
compression, a typical grid has Nr — 270 cells in the radial 
direction and Nz — 432 cells in the axial direction. The num- 
ber of grids covering the disc in the vertical direction (in the 
middle of the disc) is about 200. The grid also has sufficiently 
high resolution in the r— direction. The simulation region is 
stretched in the z— direction to allow for matter flow in the 
corona. 
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3.2 The parameters in different simulation runs 

We ran various simulations given different initial configura- 
tions of the field within the disc (see Fig. [T] and Tab.[T]). We 
observed that similar MRI-driven turbulence develops in all 
cases. We chose the tapered field configuration as a base, be- 
cause it is the simplest after the vertical field configuration, 
but most of its magnetic flux is confined within the disc. 

The strength of the MRI-driven turbulence and the subse- 
quent accretion rate depends on the value of the seed poloidal 
field in the disc, Bd. We observed that if the field is strong 
> 0.01 (see Sec.|A2]for reference units), the accretion rate 
is very high, and also the long wavelength modes dominate. 
If the field is weak, Bd < 0.0005) then the accretion rate is 
very low. We chose the intermediate value Bz = ±0.002 for 
almost all simulation runs. We also performed test runs at 
higher {Bd = -0.005) and lower {Bd = 0.001) fields. The 
result also depends on the orientation of the seed field in the 
disc. Note that the field of the star at the disc-magnetosphere 
boundary is positive, and hence the plus sign for Bd means 
that the fields have the same direction. 

Our initial setup for the disc and corona differs from the 
commonly used setup for accretion to black holes, where the 
matter is initially confined inside a thick torus of constant spe- 
cific angular momentum. The torus gradually evolves into the 



disc-type structure (e.g., l Hawley|2000 " ] ) . In contrast, our disc 
rotates with Keplerian velocity and has a quasi-equilibrium 
density distribution from the beginning of simulations (see 
Sec. |A4| ). In addition, matter with a seed field can flow inward 
from the side boundary, and this allows for longer simulation 
runs (see Sec.[A5]). 

We also performed simulations for different sizes of the 
magnetospheres which are determined by the dimensionless 
parameter fl (see Tab. [T] and Sec.|6|. 



3.3 Initial perturbations in the disc 

Here, we use our dimensionless variables (see Sec. |A2| ) to es- 
timate the number of MRI waves per thickness of the disc. The 
wavelength of most unstable, MRI-driven modes is given by 
Amri ~ 27iVA,z/^K- The full thickness of the initially isother- 
mal disc is 2h ^ 2cs/^k, where Cs — \fvj~p is the isother- 
mal sound speed in the disc. Hence, the number of waves per 
thickness of the disc is iVMRi = 2/i/Amri ~ Cs/^kva^z- In our 
simulations (in dimensionless units, see Sec. |A2| ) Cs — 0.032, 
the density in the disc is pd ~ 1, the Alfven velocity (based 
on the Bd = 0.002 field) is va,z = Bd/V^Tvp^ ^ 5.7 x 10"^. 
Substituting these values into the initial formula, we obtain 
the number of waves per disc and also the plasma parameter 
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Figure 4. Left panels: MRI-driven accretion at different moments of time for model TPSO where the seed magnetic field in the disc is parallel to 
that of the star at the disc-magnetosphere boundary. The background shows the density, lines show the poloidal field lines. Right panels: Same 
but for the case of antiparallel fields. The top row shows the flow in the whole simulation region at two moments of time t = 600 and 1000. The 
next row shows matter flow near the star for the same cases. The 3rd and 4th rows show the temporal evolution of matter flux onto the star M 
and torque on the star associated with matter, Nm, and the magnetic field, Nf for corresponding models. 
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This value of [3 corresponds approximately to the initial value 
of /3 in the middle of the disc in our simulations. Below, we 
show the development of the MRI-driven turbulence for one 
of our models. 



3.4 Development of MRI-driven turbulence in 
simulations 

Here, we take one of our models (TPSO) and show the devel- 
opment of the MRI-driven flow We take the case of the largest 
magnetosphere (in our set) to investigate whether the large 



magnetosphere can hinder the development of turbulence or 
subsequent inward accretion towards the star. 

When we imposed small random fluctuations of the an- 
gular velocity in the disc, the fluctuations grew Fig. [s] (top 
panels) shows the distribution of density p and [3 at an early 
time (t = 1). One can see that the density increases towards 
larger radii. In our initial conditions (see Sec. |A4| ), the den- 
sity may either increase or decrease towards larger radii, and 
it is regulated by a coefficient in the initial, almost Keplerian, 
angular velocity distribution: Q{r) = (1 ± 0.02)Qk, where 
the '— ' sign means the density increases outward. We chose 
this case for all our simulation runs because it corresponds 
to softer start-up conditions and the simulations last longer. 
Note that the initial distribution of the plasma parameter /3 
is not homogeneous: it varies between /3 = 20 at the inner 
edge of the disc and /3 = 2 x 10^ at the outer edge. In the 
middle of the disc, at r = 22, we have /3 = 3 x 10^ (which 
roughly corresponds to the value derived in the previous sec- 
tion) . This strong initial variation in /3 is associated with two 
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factors: (1) the density increases towards larger radii; and 
(2) the initial magnetic field in the disc is a superposition of 
the strong dipole field that dominates at the inner regions of 
the disc, and the almost homogeneous seed initial field in the 
disc. The superposition of fields gives the field B = 0.01 at the 
inner edge of the disc, B — 0.005 in the middle (at r — 22) 
and B — 0.002 at the outer edge. Note that in the 'torus' ini- 
tial conditions (e.g., Haw ley|[2000| ) the /3 also varies across 
the torus, but not as much. 

Later, we observed that the growth rate of the instabili- 
ties is on the time scale of Keplerian rotations. Fig. [s] (second 
row) shows the density and 13 distributions at time t — 100, 
which corresponds to one Keplerian rotation at r = 21 (in 
the middle of the simulation region) . One can see that chan- 
nel modes developed at this and smaller radii. The number of 
modes per thickness of the disc expected from theory (see Eq. 
|5j is A^MRi ~ 5 for the inner disc, A^mri ~ 8 - 9 for the middle 
of the disc and A^mri = 17 for the outer disc. One can see that 
at this time, the instability already developed at the inner disc 
with few modes visible in the plot, and the correct number of 
modes have developed in the middle of the disc. The external 
regions show stretching of the initial perturbations. 

The next row in Fig.jsjshows a later time, t — 300, which 
corresponds to one Keplerian rotation at the outer edge of the 
disc r 44. One can see that the channel modes have devel- 
oped in the whole disc including the outer boundary. Again, 
the number of modes roughly corresponds to the predicted 
number A^mri ~ 17 for this value of the field. Note that re- 
cently developed channel modes in the outer part of the disc 
are very ordered, while those in the inner and middle parts of 
the disc start to be mixed. 

At later times, t > 500, the channel modes mix as a re- 
sult of interacting with one another, as well as with turbu- 
lent matter that is not in channel streams. Recently, [Latter et| 
[aL] ( |2009| ) suggested that the 'turbulent mixing' is probably 
the main process which destroys the channel modes (and not 
the parasitic instabilities as suggested earlier, e.g. [Goodman] 
[& Xu[[1994l ). The bottom panels corresponding to t = 1000 
and t — 2000 show that the flow in most of the disc is very 
well mixed and is turbulent. The isolated channel streams ap- 
pear in different parts of the disc and propagate inward, but 
they do not determine the character of matter flow at the in- 
ner parts of the disc. Formation of isolated channels which 
transfer energy to the turbulent matter has been discussed 
by [L atter et al., ( ,2009 ). Most importantly, the turbulent flow 
dominates in the inner parts of the disc which are important 
for investigating the disc-magnetosphere interaction. 

The left panels of Fig. [s] show that the density is grad- 
ually redistributed in such a way that the highest density is 
in the inner parts of the disc, which corresponds to a quasi- 
stationary disc equilibrium. 
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Figure 5. Top panels: Radial distribution of matter and magnetic 
pressure, {Pyn),{Pf), and stresses, {Tm),{Tf) which were inte- 
grated across the disc, at time t = 1000. Bottom panels: Same for 
a— parameters associated with matter and magnetic stresses, (Xm ,c6f. 
Left and right panels correspond to models TP30 (parallel fields) and 
TA30 (antiparallel fields), respectively. 
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Figure 6. Top panels: Temporal evolution of matter and magnetic 
pressure, (Pm), {Pf), and stresses , (Tm), {Tf) taken at the radius 
r = 20 and integrated across the disc. Bottom panels: Temporal vari- 
ation of a— parameters associated with matter and magnetic stresses, 
am,oif. Left and right panels correspond to models TP30 (parallel 
fields) and TA30 (antiparallel fields). 



the field of the same polarity. In the case of antiparallel fields, 
the field at the disc-magnetosphere boundary is weaker ini- 
tially, and subsequently, disc matter brings in a field of the 
opposite polarity, and this influences matter flow in the in- 
ner disc and at the boundary. Below, we consider two specific 
models with parallel (TP30) and antiparallel (TA30) fields 
and discuss the results in detail. 



4 MRI-DRIVEN ACCRETION IN CASES OF PARALLEL 
AND ANTIPARALLEL FIELDS 

Our simulations show that the disc-magnetosphere interac- 
tion depends on the orientation of the seed poloidal field 
in the disc. If at the disc-magnetosphere boundary the 
B^— component of the disc field points in the same direction 
as the field of the star, then the field is enhanced at the bound- 
ary initially. Subsequently, the disc matter continues to bring 



4.1 Matter flow, matter flux and torque 

We observed that in both models, MRI-driven turbulence de- 
velops in the disc and matter flows inward. The top rows of 
Fig. [4| show accretion on large and small scales at two times 
for model TP30 with parallel fields (left side of the plot) and 
model TA30 with antiparallel fields (right side). We observed 
that in both cases, the disc moves inward due to MRI-driven 
turbulence, but it is stopped by the magnetic pressure of the 
star. Then matter is lifted above the closed magnetosphere 
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Figure 7. Top panel: Distribution of radial velocity Vr, specific matter 
flux pvr, specific magnetic torque Nf and the density p across the 
disc at radius r = 20 and time t = 1000 in model TP30. The density 
is decreased by a factor of 50 for scaling, and Nf is increased by a 
factor of 4. Bottom panel: same as above, but for three components 
of the magnetic field and plasma parameter /3 = STrPm/B'^. The 
plasma parameter varies between small values outside the disc, up to 
/3 » 1 inside the disc. 



and accretes to the stellar surface in funnel streams. We ob- 
served that the closed magnetosphere initially expands to- 
wards one side, while matter accretes from the other side. Ini- 
tial penetration of the disc matter through the closed magne- 
tosphere is accompanied by strong reconnection events. Note 
that the disc is thicker in the case of antiparallel fields. 

We calculated matter fluxes and torque at the surface of 
the star: 



dS • (Nn, + Nf ) , 



(7) 



(8) 



where dS is the surface area element directed outward; Nm 
and Nf are torques associated with matter and the magnetic 
field: 

= rpv^^rp , Nf = -r^^ . (9) 

Fig. [4| (bottom panels) shows that the matter fluxes are 
strikingly different in cases of parallel (left panels) and an- 
tiparallel (right panels) fields. In the case of parallel fields, 
matter accretes onto the star in bursts, whereas in the case of 
antiparallel fields the accretion is smooth (after the first burst) 
and is also at the lower level. We suggest that in the case 
of parallel fields, the magnetic flux accumulates at the disc- 
magnetosphere boundary and blocks accretion. Matter accu- 
mulates, and later finds a path towards the star. This leads 
to periods of low accretion rate (when accretion is blocked) 
and high accretion rate through a burst. For antiparallel fields 
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Figure 8. Top panel: close view of the funnel stream at T = 1000 in 
model TP30. Bottom panel: radial distribution of density (p), angular 
velocity (Q), plasma parameters (/3, /3i), azimuthal component of 
the field, B^, and radial and z-components of the field (Br, Bz). The 
density is multiplied by a factor of 0.1 for scaling. 



fields of opposite polarity frequently reconnect at the disc- 
magnetosphere boundary, opening the path for smoother ac- 
cretion. It is interesting to note that the 'first burst' is observed 
in both cases. 

We calculated the torque exerted on the star by matter, 
Nm and the field Nf. We observed that the magnetic torque 
is much larger than the material torque and it repeats the pat- 
tern of the matter flux, M. The torque is positive (we changed 
the sign compared with the formulae Eq.|9]for convenience of 
illustration) so the accreting matter spins up the star. This is 
what we expect, because the star rotates slowly. The observed 
properties of the magnetic and material torques are similar to 
those observed in laminar, a— discs ( |Romanova et al.,2002] ). 



4.2 Analysis of stresses and a— parameters 

Here we analyze the stresses acting inside the disc. We sep- 
arate the disc from the low-density corona using one of the 
density levels, pdisc — 0.3, which is typical for the boundary 
between the disc and corona. We integrate matter (subscript 
'm') and magnetic (subscript 'f') stresses in the z— direction 
and obtain their distribution as a function of radius in the 
co-moving frame: 
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Figure 9. Top Panels: Matter flow before, during, and after the ac- 
cretion event in model TP30, where the accretion rate enhancement 
is determined by the 'push' mechanism. The background shows the 
density and the lines are selected magnetic field lines. Bottom Panel: 
Matter flux onto the star during this accretion event, which is marked 
by an arrow. 
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Figure 10. The figure shows an event of the 'first burst', where the 
initial stage of the disc penetration leads to reconnection of the sig- 
nificant magnetic flux of the external magnetosphere. An example is 
shown for model TP20. 

^ 2h J 47T 

where 

are averaged azimuthal velocity and matter flux, and h = 
h{r) is the half- thickness of the disc. 

The matter and magnetic pressure are 

The standard a— parameters associated with matter and mag- 
netic stresses are: 

_2(7^ _ 2 (Tf) 

Fig. [s] shows the radial distribution of the magnetic and 
matter stresses, pressure, and corresponding a— parameters 
for models TP30 and TA30. One can see that in model TP30 
(see left panels) the magnetic stress (T/) is about 10 times 
larger than matter stress and hence the magnetic stress de- 
termines the angular momentum transport. The matter pres- 
sure dominates over the magnetic pressure. The magnetic 



a— parameter varies in the range of af ^ 0.01 — 0.05, while 
the matter a— parameter, am, is 5 — 10 times smaller. Hence, 
the angular momentum transport and inward accretion is 
determined by the magnetic stress. Note that the size of 
the stellar magnetosphere (magnetically-dominated region) 
is rm ~ 4 - 6i?*, and this region should be excluded from our 
consideration, because we are only interested in the stresses 
in the disc. 

The right panels of Fig. [s] show corresponding values for 
model TA30. In this model the seed poloidal field in the disc 
is antiparallel to the stellar field, and initial stress is some- 
what smaller in the inner parts of the disc. We should note 
that in this case, the field in the disc is weaker than in the 
model TP30, because the oppositely-directed stellar field is 
subtracted from the disc field. This leads to smaller magnetic 
stress, which is still a few times larger than the matter stress. 
The magnetic pressure is smaller than the matter pressure in 
the majority of the disc. The magnetic alpha-parameter af 
is about 1.5 times larger than am and the magnetic field is 
responsible for the angular momentum transport. 

In addition, we investigated variation of stresses and 
other variables in time. The averaged values are taken 
in the middle of the disc (at a radius r = 20). Fig. [6] (top 
panels) shows that the magnetic a— parameter (a/) is larger 
than the am -parameter during the whole time, and hence the 
magnetic stress determines accretion. Both a-parameters de- 
crease with time. The value of am has a similar pattern, but 
at smaller values. The bottom panels of Fig. |6] show that all 
stresses and pressure values are quite steady with time on 
average, and only slightly decrease with time. Note that af 
and am decrease with time, but this is partially because the 
stresses gradually decrease with time and partially because 
the matter pressure (Pm) increases. 

Results are shown up to t = 2000 periods of rotation at 
the inner boundary (r = 1), which corresponds to 20 periods 
of rotation at r = 22. 



5 THE DISC-MAGNETOSPHERE INTERACTION 

Here, we investigate in greater detail processes at the disc- 
magnetosphere boundary. 



5.1 Where the disc stops 

We take as an example the model TP30 and one of the mo- 
ments in time when matter flows along the funnel towards 
the star. Fig. [s] (top panel) shows the inner part of the simula- 
tion region. One can see that the disc is stopped at some dis- 
tance from the star and matter flows into the funnel stream. 
The conductivity is high, and some matter is dragged into the 
funnel stream. The matter energy-density is still higher than 
the magnetic energy in the funnel, but the magnetic field im- 
poses some tension force, which acts in the direction opposite 
to the flow. This stress is released due to reconnection inside 
the funnel stream. Such mini-reconnection events are always 
observed in the funnel streams, in particular in cases of par- 
allel fields. 

The bottom panel of Fig. [s] shows a linear distribution of 
different variables in the equatorial plane. One can see that 
at the disc-magnetosphere boundary the density p drops from 
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Figure 11. Sample plots for accretion onto stars with different magnetosphere sizes at time t = 700. Parameter jl varies (from left to right) as 
jl = 2,5,10,20,30. 



high values in the disc to very low values in the magneto- 
sphere. The component of the field dominates inside the 
disc, but drops at the boundary. The Bz and Br components 
are very small in the disc, but strongly increase in the region 
of the magnetically-dominated magnetosphere. Note that the 
Br— component is quite large because the magnetosphere is 
pushed upward by the funnel stream. Fig. [s] also shows the 
equatorial distribution of plasma parameter /3 = SnPm/B'^ 
and the modified parameter [3i based on the total material 
stress: 



Pi 



Pm + PV'^ 



(10) 



One can see that the disc-magnetosphere boundary is located 
between the radii where ^ — I and — 1. This is in accord 
with earlier simulations of a— discs, where the inner bound- 
ary usually coincides with Pi — I surface (e.g., |Romanova et] 
|aL]|2002 , Romanova et al.j 2004| ), while the funnel starts at 
13 = 1 (Bessolaz e t al.,2008X 

However, the disc-magnetosphere interaction in the case 
of turbulent discs shows many features which are different 
from those observed in the case of a— discs. One of the differ- 
ences is that in the modeUng of accretion from a— discs the 
diffusivity was taken to be of the same order as viscosity, and 
both can be small, a ^0.01 — 0.02 (like, e.g., in Romanova et 
al. 2002, 2003), or larger, a = 0.1 (e.g. Bessolaz et al.|2008l ), 
and in both cases matter smoothly flows from the disc to the 
star. The MRI-driven turbulence provides viscosity at the level 
of a ^ 0.01 — 0.1, but it does not provide the cross-field diffu- 
sivity. This is why we often have the situation that accretion 
is blocked by the magnetic flux, and it leads to strong vari- 
ability. In addition, in cases of one or another polarity in the 
disc, the variability is different, as discussed in Sec. |4.1| and 
subsequent sections. 



5.2 The 'push' mechanism of accretion and oscillations 

We often see strong oscillations in the light-curves. Usually, 
oscillations are connected with an accumulation of magnetic 
flux at the boundary, which generally occurs in cases where 
the seed magnetic field has the same polarity as the magnetic 
field of the star (cases of parallel fields) . 

Fig. |9] shows one of the episodes of enhanced accretion 
observed in model TP30 (the case of parallel fields, see Fig. 
|4]) . The left panel of Fig. [9] shows that matter accumulates 
at the disc-magnetosphere boundary. The magnetic field car- 
ried by the inner disc has the same polarity as the stellar 



field, and hence they do not reconnect. Matter bends the 
magnetosphere in such a way that accretion towards the 
star through the funnel stream becomes impossible. How- 
ever, when 'enough' matter accumulates, it starts pushing the 
magnetosphere forward, so that accretion becomes possible 
(see middle panel). This happens because the gravity force 
becomes larger than the tension force associated with bend- 
ing. After downloading matter onto the surface of the star, 
the magnetosphere expands and accretion is blocked again 
by the negative inclination of the field lines (see right panel) . 
The bottom panel shows the corresponding burst in the accre- 
tion rate. We call this mechanism the 'push' mechanism. It is 
responsible for the majority of bursts seen in the light-curve of 
model TP30 (see Fig. [4] left panels) and in many other cases 
shown below. 



5.3 The 'First Burst' 

The accretion rate in the disc may vary with time. If accretion 
rate is relatively low, then the magnetosphere expands and 
the disc is at larger radii. When accretion rate increases after a 
period of low accretion rate, the inner disc moves forward and 
compresses the magnetic flux of the star. Finally, it stops at 
the point where material and matter stresses are comparable. 
However, the compressed magnetic flux blocks accretion. In 
MRI simulations the diffusivity is very low, and the disc matter 
cannot easily penetrate through all this flux. 

Initially, our simulations corresponded to such a situa- 
tion. We observed that matter accumulates in the inner disc 
(see Fig. [To] left panel), then moves towards the star, forc- 
ing all the external flux to reconnect (see middle and right 
panels) . This event leads to a burst in the accretion rate (see 
bottom panel) . In addition, one can expect strong X-ray flare 
associated with reconnection of significant magnetic flux. 

The phenomenon of the first burst corresponds to our 
initial stage of evolution. However, such an initial burst is ex- 
pected in stars with strongly varying accretion rate. Note that 
the 'first burst' is expected in cases of both parallel and an- 
tiparallel fields (see Fig. [4]) . 



6 ACCRETION ONTO STARS WITH DIFFERENT 
MAGNETOSPHERES AND ACCRETION RATES 

Above, we described the case of the largest (in our set of mod- 
els) magnetosphere, with = 30. We also calculated accre- 
tion onto stars with smaller-sized magnetospheres at magne- 
tospheric parameters jl — 20, 10, 5, 2 (see Tab.[T]). For homo- 
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Figure 12. Accretion onto a star with a large magnetosphere (models TP20 and TA20) in cases of parallel (left panels) and antiparallel (right 
panels) fields. Bottom panels show accretion rate onto the star. 




geneity of results we use the same, tapered topology and the 
same seed magnetic field in the disc, which can be parallel 
(sign '+') or antiparallel (sign '-') to the field of the star, and 
has a z— component value Bd = ±0.002. 

In two exploratory runs, we fixed the magnetosphere at 



jl = 20, but changed the magnetic field in the disc. In one 
case, we increased the disc field by a factor of 2.5 and used 
the antiparallel orientation (Bd = 0.005). In the other case, 
we decreased the field by half and took the parallel orientation 
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Figure 14. Same as in Fig.[T2]but for models TPS and TAS where fi = 5. Note the episodic boundary layer accretion observed in model TPS. 



tt Parallel (model tp2) 




2 



accretion through the funnel 




boundary layer 
regime 


, u_l_ 







500 550 600 time 650 700 

Figure 15. Top panel: Accretion onto a star with a tiny magnetosphere 
ji = 2 and parallel fields (model TP2). Left and right panels show 
accretion before and after the onset of the boundary layer regime. 
Bottom panel: Matter flux onto the star shows the transition from the 
funnel flow accretion (at the left) to the boundary layer regime of 
accretion (at the right) . 



moved from the bottom for TP30 and from the top for TA30) 
during the whole simulation run. In cases of a smaller magne- 
tosphere, the funnel flips between the top and bottom sides, 
although this happens only a few times per simulation run. 
The density is usually enhanced at the inner part of the disc, 
excluding the case of the smallest magnetosphere at /2 = 2, 
where the disc accretes directly through the equatorial plane, 
that is, in the boundary layer regime. 

In Sec.|4|we performed an analysis of stresses for models 
RP30 and TA30 and observed that the accretion rate strongly 
oscillates in the case of parallel fields, but is much smoother 
and smaller in the case of antiparallel fields. The question 
is, how general is this behavior? We realize that initially, the 
disc is threaded not only with the seed magnetic field, which 
we introduce at the level of Bd — ±0.002, but also with the 
magnetic flux of the star, and this flux is larger in cases of 
larger magnetospheres and can dominate, or can be negligi- 
bly small, in cases of smaller magnetospheres. Either way, we 
expect that the total amount of stress in the disc is different 
in cases with different fl. Hence, with changing fl, we also 
change the total amount of stress in the disc. Below we in- 
vestigate the situation and derive general trends for the value 
and variability of the accretion rate in different cases. 



of the field. Below, we analyze results obtained in all these 
cases. 

Fig. [Tl] shows magnetospheric flow in cases of the dif- 
ferent sizes of magnetospheres. For homogeneity, we chose 
the case of parallel fields, Bd — 0.002 and show data for the 
moment in time t — 700. One can see that the size of the mag- 
netosphere decreases with jl, as expected. The funnel stream 
forms from either top or bottom side of the magnetosphere. 
In the case of the largest magnetosphere (/i = 30), we ob- 
served that the funnel formed on one side and stayed there (it 



6.1 Matter flow and accretion rates in different models 

First, we show the character of matter flow and accretion 
rates obtained in different simulation runs. In figures [T2|13| 
and [14] we show results in the same style, where the left and 
right hand sets of figures correspond to the cases of parallel 
and antiparallel fields, respectively. 

Fig. [T2I shows results of simulations for models TP20 and 
TA20 (/x = 20). One can see that the MRI-driven instability 
has developed in both cases, and in the case of parallel fields 
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Figure 16. Left panels: Example of accretion and accretion rate in the case of an intermediate-sized magnetosphere and parallel magnetic fields 
(model TPIO). Right panel: Same, but for antiparallel magnetic fields in the disc. We chose moments of time corresponding to bursting and 
smooth accretion (one per each case) . 



matter accretes through episodes of matter accumulation and 
accretion, and the 'push' mechanism of accretion is observed. 
Simulations show that the accretion rate is high and strongly 
oscillating (see left bottom panel of Fig.[T2]). In the right pan- 
els we see that in the antiparallel case, the disc is thicker and 
not as dense. Matter is accumulated near the star in a thick 
layer, and accretes onto the star, but the accretion rate is lower 
than in the case of parallel fields, and no oscillations are ob- 
served, excluding the first burst discussed in Sec. |5.3| 

Fig. [13] shows a case with a smaller magnetosphere, 
/i = 10 (models TPIO and TAIO). One can see that in cases 
of parallel fields the disc near the star is thin and the den- 
sity increases towards the star. Matter accretes through funnel 
streams. The matter flux at the star has multiple oscillations 
up to t < 600, then it increases and there are fewer oscilla- 
tions. We suggest that there are fewer oscillations because the 
magnetosphere is smaller and it is easier for matter to climb 
the magnetosphere and form a quasi-steady flow. In the case 
of antiparallel fields, the density in the disc decreases towards 
the disc-magnetosphere boundary. The disc is thicker and the 
accretion rate is lower. We analyze the stresses and discuss 
possible reasons below. 

Fig.[l4|shows the case of an even smaller magnetosphere, 
/i = 5 (models TPS and TA5). One can see that in the case of 
parallel fields, matter smoothly accretes through the funnel 
stream, up to time t ^ 1400. However, afterwards the accre- 
tion becomes bursty, because matter reaches the surface of the 
star episodically and accretes through the boundary layer (see 
left bottom panel of Fig. [14]) . The episodes of accretion are 
followed by an expansion of the magnetosphere and an accu- 
mulation of matter. Then the process repeats, giving multiple 
bursts observed in the figure. In the case of antiparallel fields 



(right panels in the figure), the density in the disc decreases 
towards the star, and the geometry is similar to that observed 
in model TAIO. The accretion rate is quasi-stationary like in 
the TAIO case. 

Fig. [15] shows an example of a simulation in the case of 
the smallest magnetosphere considered in our paper, jl — 2. 
The figure shows only the case of parallel fields (model 
TP2), where we observed the quasi-steady regime of accre- 
tion through the boundary layer (top right panel) . In the case 
of antiparallel fields, accretion is similar to that in model 
TA5, but the magnetosphere is somewhat smaller. Interest- 
ingly, when the disc accretes through the boundary layer, the 
magnetic field of the star is compressed by the disc and ob- 
tains a complex geometry. 

There are a few factors which determine the differences 
in matter flow and accretion rates seen in different models: 

• Different orientation of the disc field relative to the stellar 
field (parallel and antiparallel cases) leads to different mag- 
netic topologies at the boundary and can influence matter 
flow at the disc-magnetosphere boundary; 

• In cases of same/opposite orientations of the field, the 
magnetic stress is higher/lower at the boundary; 

• In cases with lager/smaller jl, the stellar magnetic field 
contributes more/less to the total magnetic energy of the disc, 
so that the accretion rate is different. 

Often, these factors act simultaneously, and it is difficult 
to separate one from another. To investigate the dependence 
of processes on the accretion rate, we chose the magneto- 
sphere with jl — 20, and performed two additional runs. In 
one of them, we increased the field in the disc by a factor of 
2.5: Bd — —0.005. In the case of parallel fields, we expected 
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a higher accretion rate and bursty accretion. This is why we 
chose antiparallel fields, where the situation is less clear (the 
model TA20-005). In another test case, we decreased the field 
by a factor of two: Bd = 0.001. We know that in the case of 
antiparallel fields we would again observe a very low accre- 
tion rate. Thus, we chose the case of parallel fields, where 
the lower accretion rate could yield a different result (model 
TP20+001). 

Fig. [16] shows these two cases where the left and right 
panels correspond to these test models. We observed that for 
Bd = 0.001 (left panels), the accretion rate in the disc is lower 
compared with the base case (TP20) and the density in the in- 
ner parts of the disc is lower as well. We observed that accre- 
tion is bursty during about 600 rotations, but later a steady 
funnel stream is formed. The accretion rate onto a star is a 
few times lower compared with the TP20 case. Comparison 
of these two cases shows that for the persistent bursty ac- 
cretion (observed in TP20), one needs higher accretion rate, 
at which the magnetosphere is bent and traps the accreting 
matter. Hence, the parallel field configuration is favorable for 
matter trapping and bursty accretion, but it is not the only 
condition. 

Fig. [16] (right panels) shows that a stronger magnetic 
field in the disc generates stronger MRI-driven turbulence, 
and the accretion rate is much higher than in any of the above 
models. We observed that in such a situation, the inner disc 
bends the magnetosphere; the 'push' mechanism is switched 
on and we observe the persistent bursty accretion. This is not 
typical for the case of antiparallel fields in above models. We 
conclude that the bursty accretion can appear in the cases 
of both parallel and antiparallel fields if the accretion rate 
is high. We should note, however, that this is a special case 
where the accretion rate is very high, so that the initial rapid 
accretion strongly reconstructs the disc. 

In all above cases, we have different initial magnetic flux 
in the disc. To compare the total magnetic energy stored in the 
disc, we separated the disc using the density level p = 0.3, 
integrated stresses and pressure in the whole disc, and cal- 
culated these values per unit volume. Fig. [Tt] shows values 
af = 2/3{Tf)/{Pm) and ^ = {Pm)/{Pf). The right panel 
of the figure shows that the value of /3 is very high in the 
beginning, but drops rapidly because the B^— field strongly 
increases. So, at t = 1000, varies between 3 and 6 for 
our main models, and is about 1 and 12 for test models. 
These numbers are reasonable for the MRI accretion discs 
(e.g., |Hawley ||2000|) . The volume-averaged parameter a (at 
t = 1000) varies in the range of 0.01 — 0.016 for the main 
models, and is smaller (0.006) and larger (0.035) for test mod- 
els. We should note that apart from the model TP30, all other 
models have quite similar averaged parameters in the disc, so 
that the accretion disc properties are not that different in dif- 
ferent cases. Therefore, we can use different runs for analysis 
of other characteristics of the disc-magnetosphere boundary, 
such as dependence on reconnection and the local stress de- 
pression due to different orientations of the field. The test 
cases are different from the main models, but this was the 
goal - to check the difference, and we consider and analyze 
these cases separately. Note that the volume-averaged stresses 
and a/— parameters are about twice as small compared with 
those at the disc-magnetosphere boundary (see next section) . 



6.2 Analysis of cases of parallel and antiparallel fields 

Here we perform an additional analysis of all cases with the 
aim of understanding the role of the field orientation (parallel 
versus antiparallel) in different cases. 

Sketches in Fig. [19] show the expected topologies of the 
field at the disc-magnetosphere boundary in cases of parallel 
and antiparallel fields. The top panel shows that in the case of 
unidirectional (parallel) fields of the disc and the star, mag- 
netic flux of the same polarity is accumulated at the boundary 
and traps the matter. We suggest that this might be a fac- 
tor which leads to oscillating accretion. Our simulations of 
the main models (Bd — ±0.002) show that accretion rate 
varies from case to case, but we observe bursty accretion in 
the case of parallel fields. The test cases in models TA20-005 
and TP20+ 0.001 show that the orientation of the field is not 
the only reason for bursty accretion. Nevertheless, we clearly 
see the tendency for parallel fields to favor bursty accretion. 

The right panel in Fig. [T9] shows how the topology of the 
field at the magnetosphere boundary changes as the result of 
reconnection. It is interesting to note that reconnection be- 
tween the field of the star and oppositely oriented field of 
the disc leads to the field lines which begin at the star but 
then stretch out towards the disc, and often go along the disc. 
Hence, the reconnection lets matter flow towards the deeper 
field lines of the closed magnetosphere, and therefore acts as 
an efficient diffusivity. On the other hand, the distribution of 
the field lines after reconnection has a tendency to block for- 
mation of the funnel flow towards the star. We suggest that 
this factor can be important in both decreasing the accretion 
rate and making accretion smoother. The test model TA20- 
B005 shows a different evolution, but this is a somewhat spe- 
cial case. We think that the opposite polarity of the fields is 
not a sufficient condition for smooth accretion. However, at 
moderate accretion rates, we would expect less bursty accre- 
tion than in the case of parallel fields. 

Note that the situation is analogous to the Solar wind 
interaction with the Earth's magnetosphere, where one orien- 
tation of the solar wind B -field allows for reconnection and 
the other one does not (e.g., |Frey et al|2003] ). 

Another important issue is that the magnetic stress at 
the disc-magnetosphere boundary is larger/smaller in cases 
of parallel/antiparallel fields. This can lead to a local de- 
crease in the accretion rate at the boundary. We integrated 
stresses in the z— direction, and plotted them as a function 
of radius at t = 1000. Fig. [Ts] (right panels) shows the ra- 
dial distribution of magnetic stresses ZT/), where red and 
blue lines indicate parallel and antiparallel fields, respectively. 
One can see that at the inner boundary, stresses are system- 
atically (about 10 times) larger in cases of parallel fields. The 
region where stresses are different is located to the right of 
the magnetically-dominated magnetosphere, and has a typ- 
ical thickness of Arstress = 4i?o. The size of the magneto- 
sphere decreases with jl and hence this region of different 
stresses (see yellow-color region in Fig. [Ts] ) moves slightly to 
the left. The left panel shows that the a— parameter is also 
systematically larger in cases of parallel fields, although the 
result is not as clean because matter pressure (Pm) varies and 
is in the denominator of af. For example, (Pm) is higher in 
cases of parallel fields. Nevertheless, we see enhancements of 
a in models with fl — 20, 10, 5. We conclude that the orien- 
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Figure 18. Left panels: Radial distribution of the magnetic diffusivity parameter af in the equatorial plane at sample time t = 1000 for models 
with different magnetospheric sizes (parameters p.). Right panels: same but for magnetic stresses. The bottom row shows test cases with 2.5 
times larger (model TA20-005) and 2 times smaller (model TP20+001) values of the magnetic field in the disc, compared with the base field of 
Bd = 0.002. 
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Figure 19. Top panels: Sketch of the disc-magnetosphere interaction in case of parallel fields (where the disc magnetic field has the same direction 
as the stellar field). In this case, the magnetic flux of the same polarity is accumulated at the boundary, and the magnetosphere is bent by the 
disc in the equatorial plane. Bottom panels: Same but for antiparallel fields. The dashed curve in the left panel shows the suggested field topology 
after reconnection. The dashed lines in the right panel show the X — type reconnection point. 



tation of the field does change the local accretion rate in the 
region of the inner disc that is adjacent to the magnetosphere. 



Fig. 18 also gives a representative number for a/ in 
different models. One can see that this parameter varies in 
the range of a/ ^ 0.01 — 0.04 with an average value of 
af ^ 0.02 - 0.03 for all our standard models. Note that at 
r > 10, discs with parallel and antiparallel fields give similar 
af. 

The bottom plots of Fig. [Ts] show our test cases. One can 
see that for Bd = —0.005, the stress is much larger and a/ ~ 
0.1 on average, while for Bd = 0.001, the stress is smaller 
than in our main cases, and a/ ~ 0.01 — 0.02. 



7 MAGNETICALLY-DOMINATED CORONA 

Initially, we have a low-density corona threaded by the dipole 
magnetic field. The differential rotation of the foot-points of 
loops (threading the disc and the star) could lead to the in- 
flation of the field lines and to the expansion of the corona 
([ Lovelace et aL| p!99P). Both, the coronal expansion and the 
MRI instability occur with the same time-scale of Keplerian 
rotation at a given radius r. However, we observed that 
the global inflation of the coronal magnetosphere has been 
strongly disturbed by the initial penetration of the inner disc 
through the external magnetosphere. We observed that the 
initial penetration led to a strong reconnection event (see 
Fig. [Tq| ) and a large magnetic island formed and moved out- 
ward along the disc. This island pushed the coronal field 
lines threading the disc to larger radii. Subsequent events 
of reconnection led to the formation of many more (but 
smaller) islands, which formed the turbulent component of 
the magnetically-dominated corona. The field lines connect- 
ing a star with the external corona inflated and formed an 
elongated field structure propagating up and down. Subse- 



quent evolution led to the winding of the coronal field lines 
and to the formation of a magnetically-dominated corona 
with a high azimuthal component of the field. The corona ex- 
panded slowly to larger and larger distances from the disc. We 
investigated the corona using one of our typical cases (model 
TP20). Simulations show that the magnetically-dominated 
corona reached the top and bottom boundaries of the region 
at t ^ 3000Po, which corresponds to about 30 rotations at 
the grid center (in the r— direction). The inner parts of the 
corona, closer to the star, are turbulent and represent multiple 
magnetic islands. The external 2/3 of the disc are threaded 
with the ordered magnetic field expanding outward. Interest- 
ingly, no disc wind flows from the disc to the corona along 
these lines, in spite of the high inclination of these field Hnes 
towards the disc. In addition, the magnetic field is tightly 
wrapped in the corona and a high gradient of magnetic pres- 
sure is observed. In spite of that, neither centr if ugally- driven, 
nor magnetically-driven wind was formed along these field 
lines. This issue requires further investigation. 

This magnetic, slowly moving corona resembles the 
corona formed in [Miller & Stone] ( [2000| ) simulations. How- 
ever, in their simulations, the field in the corona has been 
generated in the disc and moved upward as a result of buoy- 
ancy. In our simulations, the field originates from the initial 
coronal field of the star, but is strongly wrapped by the rota- 
tion of the disc. 

Fig. [20] shows different features of the corona. The left 
panel shows that the corona has low density. Next panel (to 
the right) shows that the expanding region of the corona is 
magnetically-dominated (plasma parameter /3 << 1). The 
next plot shows that the corona rotates approximately with 
an angular velocity of the disc. Next plot shows that the az- 
imuthal field Bcf) is strong and that there is an outwardly di- 
rected gradient of magnetic pressure. The azimuthal field is 
about 10-100 times stronger than the poloidal field. The right 
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panel shows that matter in the expanding corona is cold and 
hence the corona expands due to magnetic pressure. 

This new type of magnetic structure is not a jet but in- 
stead a magnetic formation filled with the ordered (mainly az- 
imuthal) and turbulent magnetic field which slowly expands 
in the axal direction. 



8 ASTROPHYSICAL EXAMPLES 

Table |A1] shows typical reference values. Below are a few 
examples where the values obtained in simulations are con- 
verted to dimensional values for typical CTTSs, magnetic CVs 
and milHsecond pulsars. 



8.1 Application to Classical T Tauri stars 

We consider a few simulation runs relevant to /i = 20 (TP20, 
TA20, etc.) and convert dimensionless numbers presented in 
the paper to dimensional values. 

The reference values shown in Tab. |A1] are relevant to 
stars with different jl, but the magnetic field of the star is 
different for different fl. For fl — 20, we obtain an equatorial 
field: Bl"^ = jlB^ = lOOOG (and polar field B^ = fl ^ Bo = 
2000G). The z— component of the seed magnetic field in the 
disc is Bd = 0. 002^0 = O.IG, and the azimuthal component 
of the field amplifies up to 0.5 — IG. The magnetic field 
of the star at the disc-magnetosphere boundary (say at 3R^) 
is B = 37G, that is, much larger than the field in the disc. Fig. 



12 (bottom panels) shows accretion rate M in cases of paral- 
lel and antiparallel fields. In the case of parallel fields we see 

that the accretion rate varies in the range M ^ 0.1 — 0.8 with 
bursts up to 1.5, which corresponds to dimensional values 

Mparaiiei = MMq ^ (0.6 - 3.6) X lO~^M0yr"^ with bursts 
up to 8.4 X lO"^M0yr~\ In model TA20 (with antiparallel 

fields), the accretion rate is much smaller, M = 0.01 — 0.05, 

and the dimensional accretion rate is Mantiparaiiei = MMq ^ 
(0.1 - 2.8) X lO~^M0yr~\ Hence, if the global magnetic 
field alternates polarities, then the accretion rate may vary 
from higher values and oscillating character to lower rates 
with no oscillations. The time-scale between oscillations is 
At X Po = 0.37days x 35 = 13 days, and about half of this 
value for later times. Note that values of accretion rate will in- 
crease if we take a larger magnetic field of the star and hence 
larger Bq. The dependence is M ^ B'^. 

8.2 Application to millisecond pulsars 



For /i 

Pr = 



= 20, we obtain the equatorial field of the star 
flBo = lO^G and the polar field B^ = 2 x 10^ 
G. The z— component of the seed magnetic field in the disc, 
Bd = 0. 002^0 = lO^G, and the azimuthal component of the 
field amplifies up to ^ lO^G. Fig. 12 (bottom panels) 

shows accretion rate M in cases of parallel and antiparallel 
fields. In the case of parallel fields the accretion rate varies in 

the range M ^ 0.1 — 0.8, with bursts up to 1.5, which cor- 
responds to dimensional values Mparaiiei — MMq ^ (0.23 — 
1.8) X lO~^M0yr"\ with bursts up to 3.4 x lO~^M0yr"\ 
In model TA20 (with antiparallel fields), the accretion rate is 



much smaller, M — 0.01 — 0.05, and the dimensional accre- 
tion rate is Mantipara^e^ = MMo ^ (0.23-1) xlO"^M0yr-\ 
Hence, if the global magnetic field alternates polarities, then 
the accretion rate may vary from higher values and oscillating 
character to lower rates with no oscillations. The time scale 
between oscillations is At x Pq — 0.46ms x 35 = 16 ms, 
and about half of this value for later times. Note that the in- 
terval between bursts increases with the spin of the star and 
reaches the value of 300 in the propeller regime fUstyugova et| 
al. 2011). This gives a much longer interval between bursts, 
Atbursts — 138 ms. In a real astrophysical situation, where 
the diffusivity is probably even smaller than in our simula- 
tions, one can expect even longer intervals between bursts. 
Some millisecond pulsars show accretion in a flaring regime. 
For example, the pulsar SAX J1808. 4-3658 shows flaring ac- 
tivity in the tail of the burst with a quasi-period of IHz ( |Pa-| 
[truno et aL]|2009| ) . We suggest that these flares may be con- 
nected with the disc-magnetosphere interaction in the pro- 
peller regime (see also Romanova et al.||2004| |2005[ |Ustyu-| 



Igova et al.|2006[ [Romanova et al.|2009 ) 



9 CONCLUSIONS 

We performed axisymmetric simulations of accretion onto 
magnetized stars from MRI-driven discs. The main conclu- 
sions are the following: 

• Long-lasting, MRI-driven accretion has been developed 
in the disc where the magnetic stresses determine the accre- 
tion rate with equivalent parameter am ~ 0.02 — 0.04 (with 
larger and smaller values in test models) . 

• We observed that the matter-dominated turbulent disc is 
stopped by the dipole field of the star at the radius where 
the magnetic stress of the magnetosphere matches the matter 
stress in the disc. Matter is lifted above the magnetosphere 
and flows in funnel stream, which usually picks one of the 
sides (top or bottom), and may switch sides a few times dur- 
ing the simulation run. 

• Processes at the disc-magnetosphere boundary depend 
on the orientation of the poloidal field in the disc: 

(1) . Parallel fields: If the field has the same direction 
as the stellar field at the boundary, then the magnetic stress 
and accretion rate are locally enhanced. However, the field 
of same polarity is accumulated at the boundary and matter 
is trapped. This matter bends the field lines of the magne- 
tosphere near the equator, and accretion is prohibited. Later, 
matter is accumulated and flows to the star through a tem- 
porary funnel stream, and the magnetosphere expands again. 
Then the process repeats. This leads to strong oscillations of 
the accretion flux on the star. Such oscillations appear as the 
result of a relatively high accretion rate, as determined by 
MRI, and low diffusivity at the disc-magnetosphere boundary. 
In a test run with lower accretion rate {am ~ 0.01), tempo- 
rary oscillations were followed by quasi-stationary accretion 
through the funnel flow. 

(2) . Antiparallel fields: If the fields of the star and the 
disc have opposite directions, then reconnection of the field 
at the disc-magnetosphere boundary can help matter pen- 
etrate towards deeper layers of the closed magnetosphere. 
At the same time, the magnetic flux cancellation lowers the 
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Figure 20. Distribution of different values in the whole simulation region in the model TP20 at time T = 3000. The lines are sample magnetic field 
lines. The color background (from left to right): density (p), /5— parameter, angular velocity (Q), azimuthal magnetic field (B^), and temperature 
(T). 



magnetic stress locally, within a few stellar radii of the mag- 
netosphere. This leads to lower accretion rate. In addition, 
the magnetic configuration of the stellar field after reconnec- 
tion is not favorable for formation of the funnel stream. Both 
factors lead to quasi-steady accretion through a weak funnel 
stream, and to lower accretion rates compared with the case 
of parallel fields. In the test run with much higher accretion 
rate (am ~ 0.1), the effect of stress diminishing and other 
factors becomes unimportant, and matter accretes in bursty 
regime like in case (1). 

• MRI-driven turbulence provides outward angular trans- 
port and inward flow of matter, and gives a significant viscos- 
ity. However, it does not give magnetic diffusivity at the disc- 
magnetosphere boundary. We often have a situation where 
the Prandtl number is high, Prm = i^vis/i^dif » 1, at the 
disc-magnetosphere boundary. In many cases, reconnection 
helps matter penetrate towards deeper layers of the magne- 
tosphere and acts as an efficient diffusivity. In other cases, 
small, but finite numerical diffusivity helps. 

• We point out the importance of the initial interaction of 
the disc with the magnetosphere, where the disc comes from 
some distance after a period of low accretion. Such a situation 
is expected in many astrophysical systems where a star has 
a dynamically-important magnetic field. Then the inner disc 
compresses the star's magnetic flux and moves towards the 
star. Matter, after accumulation, flows towards the star and 
pushes this magnetic flux to reconnect. This leads to a strong 
burst in the accretion rate. This may also lead to a strong X-ray 
flare resulting from the reconnection of the magnetic field. 

• A magnetically-dominated corona forms above and be- 
low the disc and slowly expands in the axial direction. The 
magnetic field is tightly wrapped and the azimuthal field is 
10-100 times stronger than the poloidal field. The corona is 
cold and it expands due to the magnetic pressure gradient. 

The simulations presented here are restricted to axisym- 
metric conditions. Earlier comparisons of axisymmetric and 
3D simulations of MRI-driven accretion discs show good sim- 



ilarity ( [Hawley|[2000| ), and we suggest that our disc is suf- 
ficient for investigation of the disc-magnetosphere interac- 
tion. However, the main restriction is probably in the fact that 
in the full 3D approach matter can penetrate the magnetic 
field due to the Raleigh-Taylor instability ( A rons & Lea|1976| ). 
Global 3D simulations show that for some conditions, the disc 
matter can penetrate the magnetosphere of the star, in par- 
ticular in cases of small tilts of the dipole field and slowly 
rotating stars ( [Romanova et al.||2005[ |Kulkarni & Romanovaj 
|2008| ) . Hence, it is important to investigate the MRI-driven ac- 
cretion onto magnetized stars in full 3D approach. We discuss 
results of global 3D simulations in [Romanova et al.] ( [2011| ). 

The described research may have multiple applications 
to different astrophysical objects. Many X-ray binaries show 
accretion in a flaring regime. For example, the pulsar SAX 
Jl 808.4-3658 shows flaring activity in the tail of the burst 
with a quasi-period of IHz ( Patruno et al.|2009) . We suggest 
that these flares may be connected with the process of mat- 
ter accumulation at the disc-magnetosphere boundary with 
subsequent accretion, similar to that observed in our case 
of parallel fields (see Fig. [9| . In this case the unidirectional 
magnetic flux blocks the accretion. This effect is similar to 
the propeller effect, where the centrifugal barrier blocks the 



accretion, and accretion occurs episodically (e.g., Romanova 
2004| |2005[ [U styugova et al. 2006! [Romanova et al 



et al. 



2009 



see also [Spruit & Taam ( 1990 , 19931 andjD'Angelo & 



Spruit[ ( [2bl0| )). In this paper we considered only slowly rotat- 
ing stars. Prehminary axisymmetric simulations of accretion 
in the propeller regime from MRI-driven discs show that the 
interval between accretion events strongly increases with the 
spin of the star ( [Ustyugova et al.|2011| ) . 
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APPENDIX A: CODE DESCRIPTION, INITIAL AND 
BOUNDARY CONDITIONS AND REFERENCE UNITS 

Al Basic Equations 

The matter flow is described by the equations of ideal MHD: 
|J + V • (pv) = , (Al) 

%^+V.r=pg, (A2) 





CTTSs 


White dwarfs 


Neutron stars 




0.8 


1 


1.4 






5000 km 


10 km 


Bi^ (G) 


5 X 10^ 


5 X 10^ 


5 X 10^ 


Ro (cm) 


1.4 X IQii 


1.4 X 10^ 


2.9 X 10^ 


vq (cm s~^) 


2.8 X 10'^ 


5.2 X 10^ 


1.4 X 10^0 


^^0 (s-i) 


2 X 10-4 


1.03 


1.4 X 10^ 


Po 


0.37 days 


6.1 s 


0.46 ms 


Bo (G) 


50 


5.0 X 10^ 


5 X 10'^ 


PO (g cm-3) 


3.3 X 10-12 


3.7 X 10-s 


5.3 X 10-5 


PO (dy cm-2) 


2.5 X 10^ 


1.0 X 10^0 


1.0 X 10^^ 


Mo(M0yr-i) 


5.6 X 10-s 


1.5 X 10-'^ 


2.3 X 10-s 


A^o (g cm^s~^ 


)6.9 X 10^6 


1.2 X 10^6 


1.0 X 10^4 


Eo (ergs-i) 


1.3 X 10^3 


1.3 X 10^6 


1.4 X 10^8 



Table Al. Sample reference units for typical CTTSs, white dwarfs, 
and neutron stars. The field B^ corresponds to the case jl = 10. Real 
dimensional values for variables can be obtained by multiplying the 
dimensionless values of variables by these reference units. 



dB 



■ V X (v X B) = 



(A3) 



+ V • {pSw) = , 



(A4) 



d{pS) 

dt 

Here, p is the density and S is the specific entropy; v is the 
flow velocity; B is the magnetic field; T is the momentum 
flux-density tensor; and g is the gravitational acceleration due 
to the star, which has mass M. The total mass of the disc is 
negligible compared to M. The plasma is considered to be an 
ideal gas with adiabatic index 7 = 5/3, and S — ln(p/p^). 
We use cylindrical coordinates (r, 2;, 0). The condition for ax- 
isymmetry is d/dcj) — 0. 



A2 Reference Units 

The MHD equations are solved using dimensionless variables 
A. To obtain the physical dimensional values A, the dimen- 
sionless values A should be multiplied by the corresponding 
reference units Aq] A — AAq, To choose the reference units, 
we first choose the stellar mass and radius R^. The ref- 
erence units are then chosen as follows: mass Mo = M^, 
distance Rq = R^, velocity vq = {GMq/RqY^'^, time scale 
Po = 27iRo/vo, angular velocity Qq = vq /Rq. The equatorial 
magnetic field is determined as — pBo where Bo is the ref- 
erence magnetic field at Ro . The parameter p is used to vary 
the magnetic field of the star. We chose Bo such a way, that 
for p = 10, we have a field typical for a given type of stars. For 
example, for CTTS, the choice of Bo = 50G, gives equatorial 
field B^ ^ p ^ 500G (and polar field lOOOG) which is typi- 
cal for CTTSs. Then the models with larger/smaller values of 
p, correspond to smaller/larger magnetic field of the star. We 
then define the reference dipole moment po — BoRo, den- 
sity Po = Bq/vI, pressure po = povi, mass accretion rate 
Mo = povoRo, the torque A^o = po'^^o^o- energy per unit time 
^0 = po'^o^o- Temperature To Tlpo/po, where K is the gas 
constant. , and the effective blackbody temperature 

Therefore, the dimensionless variables are f = r/Ro, 
V — v/vo, t — t/Po, B — B/Bo, p — p/ Po and so on. In 
the subsequent sections, we show dimensionless values for 



all quantities and drop the tildes (^) . Our dimensionless sim- 
ulations are applicable to different astro physical objects with 
different scales. We list the reference values for typical CTTSs, 
cataclysmic variables, and millisecond pulsars in Tab.|Al[ 



A3 Code description 

We solve the full set of axisymmetric ideal MHD equations 
in cylindrical coordinates using the Godunov-type numerical 
method. To solve the set of ideal MHD equations, we use an 
approximate multi-state HLL Riemann solver similar to one 
descri bed recently by |Miyoshi & Kusano| ( [2005| ) . Compared 
with Miyoshi & Kusano ( 2005 1), we solve an equation for the 
entropy instead of full energy equation. This approximation 
is valid in cases where shocks are not important, like in our 
case. In Miyos hi & Kusano| ( |2005 ), the procedure of the cal- 
culations of the MHD equations is described in detail. Here, 
we briefly summarize the main features of our method: the 
MHD -variables are calculated in four states bounded by five 
MHD discontinuities: the contact discontinuity, two Alfven 
waves and two fast magnetosponic waves. Compared with 
[Miyoshi & Kusano] ( [2005| ), we performed the procedure of 
correction for the fast magnetosonic waves velocity so that to 
support the gap between these waves and the Alfven waves 
which propagate behind the fast magnetosonic waves. This 
helps us to escape the "non-physical" solutions of the Rie- 
mann problem. To satisfy the condition of the zero divergency 
of the magnetic field, we introduced the 0— component of the 
magnetic field potential which has been calculated using the 



method proposed by [Gardiner & Stone| ( [2005| ) . To decrease 
the error in the approximation of the Lorentz force, we use 
the splitting of the field to the dipole and calculated compo- 
nents, B — Bdip + B\ and dropped terms of the order of 
B'^^p which give zero input into the Maxwellian stress ten- 
sor ( |Tanaka|1994| ) . No viscosity or diffusivity has been added 
to the code, and hence we investigated accretion driven only 
by the MRI-turbulence. Our code has been extensively tested 
(see tests and other details in |Koldoba et al.|20TlD . 
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A4 Initial conditions 



The initial conditions are similar to those used in iRomanoval 
|et al.| ( [2005] ) and [Ustyugova et aL| ( |2006| ). Here, we summa- 
rize these conditions. 

The disc is cold with temperature Td and dense with 
density pd, while the corona is hot Tc = lOOOT^ and rar- 
efied pc = O.OOlpd. We place both the disc and the corona 
into the simulation region. We assume that the initial flow is 
barotropic with p = p{p), and that there is no pressure jump 
at the boundary between the disc and corona. Then the initial 
density distribution (in dimensionless units) is the following: 



p/lZTd , p > Pb and r > n ■ 
p/lZTc , p <Ph or r <rh 



where ph is the pressure on the surface which separates the 
cold matter of the disc from the hot matter of the corona. 
On this surface the density jumps from ph/Td to ph/Tc. Here 
Th is the inner disc radius. Because the density distribution is 
barotropic, the angular velocity is constant on coaxial cylin- 
drical surfaces about the z— axis. Consequently the pressure 
distribution may be determined from the Bernoulli equation, 



F{p) + $ + $c 



const . 



Here, $ = —GM/\y\ is gravitational potential, $c = 
/r^ine ^^(0^^^ is centrifugal potential, which depends only 
on the cylindrical radius r sin 6, and 



F{v) 



1ZTd\n{p/pb) , p > Pb and rsin^ > rb 
7lTc\n{p/pb) , p <Pb or rs'mO<rb. 



The angular velocity of the disc is slightly super-Keplerian, 
Q(^z = 0) = kQk = 1 + 0.02), due to which the density and 
pressure increase towards the periphery. Inside the cylinder 
r < n the matter rotates rigidly with angular velocity of the 
star Q(rb) = ^{GM/riy^'^. We place the inner edge of the 
disc at rb = 10 and rotate a star with angular velocity Q{rb) = 
10~^/^ ^ 0.03. The disc evolves, but we keep a star rotating 
slowly to be sure that it represents the gravitational well for 
the incoming matter (see the case of rapidly rotating stars in 
[Ustyugova et al.|201l| ). 

Initially, the region is threaded with the dipole magnetic 
field of the star. We use different values of the magnetospheric 
parameter jl from 30 to 2 (see Tab.[T]). This parameter is re- 
sponsible for the size of the magnetosphere in dimensionless 
runs, in the disc. The initial density in the disc varies from 
Pd = 1 (at the inner edge) up to />cz = 4.6 (at the outer bound- 
ary). 



A5 Boundary Conditions 

On the star: The dipole magnetic field is frozen into the sur- 
face of the star, that is the normal to the surface component 
of the field is fixed, while azimuthal component has "free" 
boundary condition, dB(p/dn — 0. There is also free bound- 
ary conditions for all other variables A, dA/dn — 0. Addi- 
tional condition is applied to the poloidal velocity which for- 
bids matter to flow out of the star. In addition, in the last grid 
we adjust the total velocity vector to be parallel to the total 
magnetic field vector. This helps to strengthen the frozen-in 
condition near the star and to enhance the divergency-free 



condition at the boundary. The top and bottom boundaries: 
The values of density p, entropy s and azimuthal field B^p are 
fixed. Free conditions for all other variables {Br, Bz, Vr, Vz, 
Vci)).ln addition, Vz = Oin case if matter flows into the simula- 
tion region. The side boundary: The side boundary is divided 
into the "the disc region" at \z\ < Zd and "corona region" at 
\z\ > Zd, where 



Zd 



GM 

Fc {Rout 



where Fc{Rout) = aGM/Rout is centrifugal potential and 
Rout is the external radius. In the disc (|z| < Zd) we have an 
inflow condition for the matter. Namely, we push matter into 
the region with small velocity: 



= -5 



2 pVK{Rout) ' 



5 = 0.02, 



and with poloidal magnetic field corresponding to the initial 
magnetic field at r = Rout- We also fix p, s and B^, while 
condition for Vz is free condition. In the corona, \z\ < Zd con- 
ditions are similar to those on top and bottom boundaries. 
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